library(xlsx)
library(plyr)
library(ggpubr)
library(ggridges)

setwd("E://OneDrive//研究//8-学位论文//数据分析处理//第6章")

#数据读取
rm(list=ls())

grow_leaf <- read.xlsx("20230216.xlsx",sheetName = "Sheet1", encoding = "UTF-8")  #生长和叶片数据
grow_leaf$T <- as.factor(grow_leaf$T)
grow_leaf$L <- as.factor(grow_leaf$L)
grow_leaf$C <- as.factor(grow_leaf$C)
grow_leaf$S <- as.factor(grow_leaf$S)

grow_leaf <- grow_leaf[grow_leaf$L %in% c(0,3,4),]

cnp <- read.xlsx("20230216.xlsx",sheetName = "CNP", encoding = "UTF-8")   #元素数据
cnp$T <- as.factor(cnp$T)
cnp$L <- as.factor(cnp$L)
cnp$C <- as.factor(cnp$C)
cnp$S <- as.factor(cnp$S)

cnp <- cnp[cnp$L %in% c(0,3,4),]

#temperature <- read.xlsx("Temperature.xlsx", sheetName = "Sheet1")

light <- read.xlsx("light.xlsx",sheetName = "Sheet2", encoding = "UTF-8")  #生长和叶片数据
light2 <- data.frame()
for (i in c(1:30)) {
  T <- rep(light[1,i+1],401)
  L <- rep(light[2,i+1],401)
  re <- rep(light[3,i+1],401)
  wave <- light$NA.[c(-1,-2,-3)]
  int <- light[c(-1,-2,-3),i+1]
  light2 <- rbind(light2, data.frame(T, L, re, wave, int))
  rm(T, L, re, wave, int)
}
light2$wave <- as.numeric(light2$wave)
light2$int <- as.numeric(light2$int)

#函数
relative.yb <- function(data, Col_DV, Col_IV, Col_group = NULL) {
  df <- data
  if (!is.null(Col_group)) {
    df <- unite(df, "group", all_of(Col_group), sep = "_", remove = F)
    df$group <- as.factor(df$group)
    group <- unique(df$group)
  } else {
    df$group <- as.factor(1)
    group <- unique(df$group)
  }
  
  data_list <- list()
  out_list <- list()
  out <- data.frame()
  for (i in c(1:length(group))) {
    data_list[[i]] <- df[df$group == group[i],]
    out_list[[i]] <- df[df$group == group[i], c(Col_IV, Col_group)]
    for (j in c(1:length(Col_IV))) {
      if(j == 1) {
        CK_row <- which(data_list[[i]][,Col_IV[j]] == 0)
      } else {
        CK_row <- intersect(CK_row, which(data_list[[i]][,Col_IV[j]] == 0))
      }
    }
    for (j in c(1:length(Col_DV))) {
      x0 <- mean(data_list[[i]][CK_row,Col_DV[j]])
      relative <- (data_list[[i]][,Col_DV[j]] - x0)/x0 *100
      out_list[[i]] <- cbind(out_list[[i]], relative)
      names(out_list[[i]])[names(out_list[[i]]) == "relative"] <- Col_DV[j]
    }
    out <- rbind(out, out_list[[i]])
  }
  return(out)
}

ANOVA.yb <- function(df, Col_DV, Col_IV, Col_group = NULL) {
  library(tidyr)
  library(car)
  library(plyr)
  library(dplyr)
  library(agricolae)
  
  if (!is.null(Col_group)) {
    df <- unite(df, "group", all_of(Col_group), sep = "_", remove = F)
    df$group <- as.factor(df$group)
    group <- unique(df$group)
  } else {
    df$group <- as.factor(1)
    group <- unique(df$group)
  }
  
  
  result <- list()
  out <- list()
  
  for (i in c(1:length(group))) {
    result[[i]] <- list()
    out[[i]] <- list()
    df1 <- df[df$group == group[i],]
    df1 <- unite(df1, "IV", all_of(Col_IV), sep = "_", remove = F)
    df1$IV <- paste("IV", df1$IV, sep = "_")
    df1$IV <- as.factor(df1$IV)
    IV <- unique(df1$IV)
    
    for (j in c(1:length(Col_DV))) {
      result[[i]][[j]] <- list()
      length2 <- function (x, na.rm=T) {
        if (na.rm) sum(!is.na(x))
        else       length(x)
      }
      datac <- ddply(df1, c(Col_IV, "IV"), .drop=T,
                     .fun = function(xx, col) {
                       c(N    = length2 (xx[[col]], na.rm=T),
                         mean = mean    (xx[[col]], na.rm=T),
                         sd   = sd      (xx[[col]], na.rm=T),
                         min  = min     (xx[[col]], na.rm=T),
                         max  = max     (xx[[col]], na.rm=T)
                       )
                     },
                     Col_DV[j]
      )
      datac$se <- datac$sd / sqrt(datac$N) 
      datac$down <- datac$mean - datac$se
      datac$up <- datac$mean + datac$se
      result[[i]][[j]][[3]] <- datac
      
      model <- paste(Col_DV[j], paste(Col_IV, collapse = "*"), sep = "~")
      normal <- data.frame(group = "", DV = "", IV = "", KS = "", KS_P = "", SW = "", SW_P = "")[-1,]
      for (k in c(1:length(IV))) {
        sample <- df1[df1$IV == IV[k], Col_DV[j]]
        ks <- ks.test(sample, pnorm, mean(sample), sd(sample))
        if (sd(sample) != 0) {
          sw <- shapiro.test(sample)
        } else {
          sw <- list()
          sw[[1]] <- 0
          sw[[2]] <- 0
        }
        normal <- rbind(normal, c(as.character(group[i]), as.character(Col_DV[j]), as.character(IV[k]), 
                                  round(ks[[1]], 4), round(ks[[2]], 4), 
                                  round(sw[[1]], 4), round(sw[[2]], 4)))
      }
      names(normal) <- c("group", "DV", "IV", "KS", "KS_P", "SW", "SW_P")
      result[[i]][[j]][[1]] <- as.data.frame(leveneTest(eval(parse(text = model)), data = df1, center = mean))
      result[[i]][[j]][[1]] <- cbind(rep(group[i], length(result[[i]][[j]][[1]][,1])),
                                     rep(Col_DV[j], length(result[[i]][[j]][[1]][,1])),
                                     result[[i]][[j]][[1]])
      names(result[[i]][[j]][[1]])[1:2] <- c("group", "DV")
      result[[i]][[j]][[2]] <- as.data.frame(summary(aov(eval(parse(text = model)), data = df1))[[1]])
      result[[i]][[j]][[2]] <- cbind(rep(group[i], length(result[[i]][[j]][[2]][,1])),
                                     rep(Col_DV[j], length(result[[i]][[j]][[2]][,1])),
                                     row.names(result[[i]][[j]][[2]]),
                                     result[[i]][[j]][[2]])
      names(result[[i]][[j]][[2]])[1:3] <- c("group", "DV", "source")
      SNK <- SNK.test(aov(eval(parse(text = paste(Col_DV[j], "IV", sep = "~"))), data = df1), "IV")$group
      SNK <- data.frame(IV = rownames(SNK), SNK = SNK$groups)
      LSD <- LSD.test(aov(eval(parse(text = paste(Col_DV[j], "IV", sep = "~"))), data = df1), "IV", p.adj="none")$group
      LSD <- data.frame(IV = rownames(LSD), LSD = LSD$groups)
      duncan <- duncan.test(aov(eval(parse(text = paste(Col_DV[j], "IV", sep = "~"))), data = df1), "IV")$group
      duncan <- data.frame(IV = rownames(duncan), duncan = duncan$groups)
      HSD <- HSD.test(aov(eval(parse(text = paste(Col_DV[j], "IV", sep = "~"))), data = df1), "IV")$group
      HSD <- data.frame(IV = rownames(HSD), HSD = HSD$groups)
      mc <- merge(merge(SNK, LSD, by = "IV"), merge(duncan, HSD, by = "IV"), by = "IV")
      result[[i]][[j]][[3]] <- merge(result[[i]][[j]][[3]], merge(normal, mc, by = "IV"), by = "IV")
      result[[i]][[j]][[3]] <- result[[i]][[j]][[3]][, c("group", "DV", Col_IV, "N", "mean", "sd", "min", "max", "se", "down", "up", "KS", "KS_P", "SW", "SW_P", "SNK", "LSD", "duncan" ,"HSD")]
      result[[i]][[j]][[3]][nchar(result[[i]][[j]][[3]]$SNK) > 2,"SNK"] <- sub('^(.).*(.)$', '\\1-\\2', result[[i]][[j]][[3]]$SNK)[nchar(result[[i]][[j]][[3]]$SNK) > 2]
      result[[i]][[j]][[3]][nchar(result[[i]][[j]][[3]]$LSD) > 2,"LSD"] <- sub('^(.).*(.)$', '\\1-\\2', result[[i]][[j]][[3]]$LSD)[nchar(result[[i]][[j]][[3]]$LSD) > 2]
      result[[i]][[j]][[3]][nchar(result[[i]][[j]][[3]]$duncan) > 2,"duncan"] <- sub('^(.).*(.)$', '\\1-\\2', result[[i]][[j]][[3]]$duncan)[nchar(result[[i]][[j]][[3]]$duncan) > 2]
      result[[i]][[j]][[3]][nchar(result[[i]][[j]][[3]]$HSD) > 2,"HSD"] <- sub('^(.).*(.)$', '\\1-\\2', result[[i]][[j]][[3]]$HSD)[nchar(result[[i]][[j]][[3]]$HSD) > 2]
    }
    for (j in c(1:length(Col_DV))) {
      if (j == 1) {
        out[[i]][[1]] <- result[[i]][[j]][[1]]
        out[[i]][[2]] <- result[[i]][[j]][[2]]
        out[[i]][[3]] <- result[[i]][[j]][[3]]
      } else {
        out[[i]][[1]] <- rbind(out[[i]][[1]], result[[i]][[j]][[1]])
        out[[i]][[2]] <- rbind(out[[i]][[2]], result[[i]][[j]][[2]])
        out[[i]][[3]] <- rbind(out[[i]][[3]], result[[i]][[j]][[3]])
      }
    }
  }
  return(out)
} #方差分析；输入(列名)：df-数据表(数据框)、DV-因变量(向量)、VI-自变量/固定因子(向量)、group-多重比较分组(向量)
#输入：result[[group]][[DV]][[1-levene(n>50-KS); 2-anova; 3-summary & normality_test & multiple_comparisons]]

analysis.yb <- function(data, Col_DV, Col_IV, Col_group, Col_species, species) {
  library(tidyr)
  df <- data
  df <- unite(df, "species", all_of(Col_species), sep = "_", remove = F)
  df$species <- factor(df$species, levels = species)
  #species <- unique(df$species)
  df1 <- list()   #df1未经转换的数据（[[1]]为全部数据，[[2]]之后为按物种分类后的相对值）
  df1[[1]] <- df[,c(Col_species, Col_group, Col_IV, Col_DV)]
  for (i in c(1:length(species))) {
    df_species <- df[df[,Col_species] == species[i], c(Col_group, Col_IV, Col_DV)]
    df1[[i+1]] <- cbind(rep(species[i], length(df_species[,1])),
                        relative.yb(data = df_species, Col_DV, Col_IV, Col_group))
    names(df1[[i+1]])[1] <- Col_species
  }
  df2 <- list()
  for (i in c(1:length(df1))) {
    if (i == 1) {
      df2[[i]] <- ANOVA.yb(df = df1[[i]], Col_DV, Col_IV = c(Col_IV, Col_group, Col_species))
    } else {
      df2[[i]] <- ANOVA.yb(df = df1[[i]], Col_DV, Col_IV = c(Col_IV, Col_group))
    }
  }
  return(df2)
} #只用于博士毕业论文，未做兼容性设计 输出1原始值方差分析，2以后不同物种的方差分析

fig.yb <- function(data, Col_IV, x1, x2, Col_DV = "DV", y1, y2, 
                   Col_group = "group", g1, g2, pointrange = T) {
  library(ggplot2)
  library(tidyr)
  
  color2 <- c("#1f78b4", "#33a02c")
  color4 <- c("#1f78b4", "#7ab5e2", "#33a02c", "#bfdb76")
  
  data <- unite(data, "Group", all_of(Col_group), sep = "_", remove = F)
  for (i in c(1:length(g1))) {
    data$Group[data$Group == g1[i]] <- g2[i]
  }
  data$Group <- factor(data$Group, levels = g2)
  
  data <- unite(data, "IV", all_of(Col_IV), sep = "_", remove = F)
  data$IV <- paste("IV", data$IV, sep = "_")
  for (i in c(1:length(x1))) {
    data$IV[data$IV == x1[i]] <- x2[i]
  }
  data$IV <- factor(data$IV, levels = x2)
  
  Fig <- list()
  
  if (pointrange) {
    for (i in c(1:length(y1))) {
      Fig[[i]] <- ggplot(data = data[data[,Col_DV] == y1[i],])+
        geom_hline(yintercept = 0, linetype = "dashed")+
        geom_pointrange(aes(x = IV, y = mean, ymin = down, ymax = up, color = Group),
                        position = position_dodge(.75), size = 1)+
        geom_text(aes(x = IV, y = up, label = duncan, color = Group), angle = 90,
                  position = position_dodge2(.75), hjust = -0.5, size = 4, show.legend = F)+
        scale_color_manual(values = color2)+
        scale_fill_manual(values = color2)+
        scale_y_continuous(labels=function(x) paste0(x,"%"),
                           limits = c(min(data[data[,Col_DV] == y1[i],]$down), 
                                      1.2*max(data[data[,Col_DV] == y1[i],]$up) - 0.2*min(data[data[,Col_DV] == y1[i],]$down)))+
        labs(y = y2[i])+
        theme_bw()+
        theme(legend.title=element_blank(),
              legend.text = element_text(size = 12, face = "bold"),
              axis.text.x = element_text(size = 12, face = "bold"),
              axis.title.x = element_blank(),
              axis.title.y = element_text(size = 14, face = "bold"),
              axis.text.y = element_text(size = 12))
    }
  } else {
    for (i in c(1:length(y1))) {
      Fig[[i]] <- ggplot(data = data[data[,Col_DV] == y1[i],])+
        geom_col(aes(x = IV, y = mean, color = Group, fill = Group), width = 0.75, position = position_dodge())+
        geom_errorbar(aes(x = IV, ymin = down, ymax = up, color = Group), 
                      width = 0.2, position = position_dodge(.75), size = 1)+
        geom_text(aes(x = IV, y = up, label = duncan, color = Group), angle = 90,
                  position = position_dodge2(.75), hjust = -0.5, size = 4, show.legend = F)+
        scale_color_manual(values = color4)+
        scale_fill_manual(values = color4)+
        scale_y_continuous(expand = c(0,0),
                           limits = c(0, max(data[data[,Col_DV] == y1[i],]$up)*1.2))+
        labs(y = y2[i])+
        theme_bw()+
        theme(legend.title=element_blank(),
              legend.text = element_text(size = 12, face = "bold"),
              axis.text.x = element_text(size = 12, face = "bold"),
              axis.title.x = element_blank(),
              axis.title.y = element_text(size = 14, face = "bold"),
              axis.text.y = element_text(size = 12))
    }
  }
  return(Fig)
}

temperature.yb <- function(data) {
  library(ggplot2)
  library(gridExtra)
  library(ggpubr)
  data$d <- data$T - data$CK
  Fig_Tsum_data <- data.frame(A = c("环境温度(°C)", "增温温度(°C)", "温差(°C)", "t值", "P值"), 
                              B = c(paste(round(mean(data$CK),2), "±", 
                                          round(sd(data$CK)/sqrt(length(data$CK)),2)),
                                    paste(round(mean(data$T),2), "±", 
                                          round(sd(data$T)/sqrt(length(data$T)),2)),
                                    paste(round(mean(data$d),2), "±", 
                                          round(sd(data$d)/sqrt(length(data$d)),2)),
                                    round(t.test(data$CK, data$T, pair=TRUE)[["statistic"]],3),
                                    round(t.test(data$CK, data$T, pair=TRUE)[["p.value"]],3)))
  if (Fig_Tsum_data[5,2] < 0.001) {
    Fig_Tsum_data[5,2] <- "<0.001"
  }
  g_data <- data.frame(date = rep(data$date,2),
                       temperature = c(data$CK, data$T),
                       group = c(rep("环境对照", length(data$date)), 
                                 rep("增温处理", length(data$date))))
  Fig_T <- ggplot(g_data)+
    geom_line(aes(x = date, y = temperature, color = group), size = 1)+
    scale_x_date(date_breaks = "1 week",
                 date_labels = "%Y/%m/%d")+
    scale_y_continuous(limits = c(min(g_data$temperature), 
                                  1.25*max(g_data$temperature) - 0.25*min(g_data$temperature)))+
    labs(x = "", y = "温度(°C)")+
    scale_color_manual(values = c("#1f78b4", "#33a02c"))+
    annotation_custom(grob = tableGrob(t(Fig_Tsum_data), rows = NULL, cols = NULL),
                      ymin = 1.1*max(g_data$temperature) - 0.1*min(g_data$temperature),
                      ymax = 1.2*max(g_data$temperature) - 0.2*min(g_data$temperature))+
    theme_bw()+
    theme(legend.position="top",
          legend.title = element_blank(),
          legend.text = element_text(size = 12, face = "bold"),
          axis.text.x = element_text(size = 12, face = "bold"),
          axis.title.x = element_blank(),
          axis.title.y = element_text(size = 14),
          axis.text.y = element_text(size = 12))
  return(Fig_T)
}

#计算的数据
grow_leaf$TB <- grow_leaf$AB + grow_leaf$UB  #总生物量
grow_leaf$RS <- grow_leaf$UB / grow_leaf$AB  #根冠比
grow_leaf$SLA <- (grow_leaf$LA / 100) / (grow_leaf$LM) #比叶面积
cnp$CN <- cnp$OC / cnp$TN * 1000
cnp$CP <- cnp$OC / cnp$TP * 1000
cnp$NP <- cnp$TN / cnp$TP

#分析
Results_grow_leaf <- analysis.yb(data = grow_leaf, 
                                 Col_DV = c("Ht", "RL", "BS", "TB", "RS", "LN", "LA", "LM", "SLA", "Fv.Fo", "Fv.Fm", "Chl", "ChlN"), 
                                 Col_IV = c("T", "L"), 
                                 Col_group = "C", Col_species = "S", species = c("I", "N"))
Results_cnp <- analysis.yb(data = cnp, 
                           Col_DV = c("CN", "CP", "NP"), 
                           Col_IV = c("T", "L"), 
                           Col_group = "C", Col_species = "S", species = c("I", "N"))

#结果输出表格
sheetname <- data.frame(c("F", "Sum"),
                        c("I-F", "I-Sum"),
                        c("N-F", "N-Sum"))
if (file.exists(paste("Results ", Sys.Date(), ".xlsx", sep = ""))) {
  file.remove(paste("Results ", Sys.Date(), ".xlsx", sep = ""))
}
for (i in 1:3) {
  for (j in 2:3) {
    if (i == 1 & j == 1) {
      write.xlsx(rbind(Results_grow_leaf[[i]][[1]][[j]], Results_cnp[[i]][[1]][[j]]),
                 file = paste("Results ", Sys.Date(), ".xlsx", sep = ""),
                 sheetName = as.character(sheetname[j-1,i]), append = FALSE, row.names = FALSE)
    } else {
      write.xlsx(rbind(Results_grow_leaf[[i]][[1]][[j]], Results_cnp[[i]][[1]][[j]]),
                 file = paste("Results ", Sys.Date(), ".xlsx", sep = ""),
                 sheetName = as.character(sheetname[j-1,i]), append = TRUE, row.names = FALSE)
    }
  }
}

#结果输出图片
Fig <- c(fig.yb(data = Results_grow_leaf[[1]][[1]][[3]], 
                Col_IV = c("T", "L"), 
                x1 = c("IV_0_0", "IV_0_3", "IV_0_4", "IV_1_0", "IV_1_3", "IV_1_4"), 
                x2 = c("T0L0", "T0L1", "T0L2", "T1L0", "T1L1", "T1L2"), 
                Col_DV = "DV", 
                y1 = c("Ht", "RL", "BS", "TB", "RS", "LN", "LA", "LM", "SLA", "Fv.Fo", "Fv.Fm", "Chl", "ChlN"), 
                y2 = c("株高(cm)", "根长(cm)", "基径(mm)", "总生物量(g)", "根冠比", 
                       "叶片数", expression("叶面积 (" * mm ^ 2 * ")"), "叶干重(g)", 
                       expression("比叶面积 (" * cm ^ 2 * "/g)"), "Fv/Fo", "Fv/Fm", "叶绿素相对含量(SPAD)", "叶氮(mg/g)"), 
                Col_group = c("C", "S"), 
                g1 = c("Mix_I","Mix_N", "Mono_I", "Mono_N"), 
                g2 = c("混合种植-鬼针草","混合种植-金盏银盘", "单一种植-鬼针草", "单一种植-金盏银盘"), 
                pointrange = F),
         fig.yb(data = Results_cnp[[1]][[1]][[3]], 
                Col_IV = c("T", "L"), 
                x1 = c("IV_0_0", "IV_0_3", "IV_0_4", "IV_1_0", "IV_1_3", "IV_1_4"), 
                x2 = c("T0L0", "T0L1", "T0L2", "T1L0", "T1L1", "T1L2"), 
                Col_DV = "DV", 
                y1 = c("CN", "CP", "NP"), 
                y2 = c("C:N", "C:P", "N:P"), 
                Col_group = c("C", "S"), 
                g1 = c("Mix_I","Mix_N", "Mono_I", "Mono_N"), 
                g2 = c("混合种植-鬼针草","混合种植-金盏银盘", "单一种植-鬼针草", "单一种植-金盏银盘"), 
                pointrange = F))

Fig_I <- c(fig.yb(data = Results_grow_leaf[[2]][[1]][[3]], 
                  Col_IV = c("T", "L"), 
                  x1 = c("IV_0_0", "IV_0_3", "IV_0_4", "IV_1_0", "IV_1_3", "IV_1_4"), 
                  x2 = c("T0L0", "T0L1", "T0L2", "T1L0", "T1L1", "T1L2"), 
                  Col_DV = "DV", 
                  y1 = c("Ht", "RL", "BS", "TB", "RS", "LN", "LA", "LM", "SLA", "Fv.Fo", "Fv.Fm", "Chl", "ChlN"), 
                  y2 = c("株高", "根长", "基径", "总生物量", "根冠比", 
                         "叶片数", "叶面积", "叶干重", 
                         "比叶面积", "Fv/Fo", "Fv/Fm", "叶绿素相对含量", "叶氮"), 
                  Col_group = c("C"), 
                  g1 = c("Mix", "Mono"), 
                  g2 = c("混合种植", "单一种植"), 
                  pointrange = T),
           fig.yb(data = Results_cnp[[2]][[1]][[3]], 
                  Col_IV = c("T", "L"), 
                  x1 = c("IV_0_0", "IV_0_3", "IV_0_4", "IV_1_0", "IV_1_3", "IV_1_4"), 
                  x2 = c("T0L0", "T0L1", "T0L2", "T1L0", "T1L1", "T1L2"), 
                  Col_DV = "DV", 
                  y1 = c("CN", "CP", "NP"), 
                  y2 = c("C:N", "C:P", "N:P"), 
                  Col_group = c("C"), 
                  g1 = c("Mix", "Mono"), 
                  g2 = c("混合种植", "单一种植"), 
                  pointrange = T))

Fig_N <- c(fig.yb(data = Results_grow_leaf[[3]][[1]][[3]], 
                  Col_IV = c("T", "L"), 
                  x1 = c("IV_0_0", "IV_0_3", "IV_0_4", "IV_1_0", "IV_1_3", "IV_1_4"), 
                  x2 = c("T0L0", "T0L1", "T0L2", "T1L0", "T1L1", "T1L2"), 
                  Col_DV = "DV", 
                  y1 = c("Ht", "RL", "BS", "TB", "RS", "LN", "LA", "LM", "SLA", "Fv.Fo", "Fv.Fm", "Chl", "ChlN"), 
                  y2 = c("株高", "根长", "基径", "总生物量", "根冠比", 
                         "叶片数", "叶面积", "叶干重", 
                         "比叶面积", "Fv/Fo", "Fv/Fm", "叶绿素相对含量", "叶氮"), 
                  Col_group = c("C"), 
                  g1 = c("Mix", "Mono"), 
                  g2 = c("混合种植", "单一种植"), 
                  pointrange = T),
           fig.yb(data = Results_cnp[[3]][[1]][[3]], 
                  Col_IV = c("T", "L"), 
                  x1 = c("IV_0_0", "IV_0_3", "IV_0_4", "IV_1_0", "IV_1_3", "IV_1_4"), 
                  x2 = c("T0L0", "T0L1", "T0L2", "T1L0", "T1L1", "T1L2"), 
                  Col_DV = "DV", 
                  y1 = c("CN", "CP", "NP"), 
                  y2 = c("C:N", "C:P", "N:P"), 
                  Col_group = c("C"), 
                  g1 = c("Mix", "Mono"), 
                  g2 = c("混合种植", "单一种植"), 
                  pointrange = T))

light_fig <- ddply(light2, .(L, wave), summarize, int = mean(int))
light_fig <- light_fig[light_fig$L %in% c(0,3,4),]
light_fig$L[light_fig$L == 0] <- "L0 (CK)"
light_fig$L[light_fig$L == 3] <- "L1 (R)"
light_fig$L[light_fig$L == 4] <- "L2 (FR)"

fig6.1 <- ggplot(light_fig)+
  geom_col(aes(x = wave, y = int, fill = wave, colour = wave), show.legend = FALSE)+
  scale_fill_gradientn(values = seq(0,1,0.05),
                       colours = c('#610061','#8300b5','#6a00ff','#0000ff','#007bff','#00d5ff','#00ff92',
                                   '#36ff00','#81ff00','#c3ff00','#ffff00','#ffbe00','#ff7700','#ff2100',
                                   '#ff0000','#ff0000','#ff0000','#db0000','#b50000','#8d0000','#610000'))+
  scale_color_gradientn(values = seq(0,1,0.05),
                        colours = c('#610061','#8300b5','#6a00ff','#0000ff','#007bff','#00d5ff','#00ff92',
                                    '#36ff00','#81ff00','#c3ff00','#ffff00','#ffbe00','#ff7700','#ff2100',
                                    '#ff0000','#ff0000','#ff0000','#db0000','#b50000','#8d0000','#610000'))+
  facet_grid(L ~ .)+
  scale_x_continuous(expand = c(0,0))+
  scale_y_continuous(expand = c(0,0), limits =c(0,390))+
  labs(x = "波长(nm)", y = expression("光照强度(μW/  " * cm ^ 2 * "/nm)"))+
  theme_bw()+
  theme(axis.title = element_text(size = 14),
        axis.text = element_text(size = 12),
        strip.text = element_text(size=12, face="bold"))
                                    
fig6.2 <- ggarrange(Fig_I[[1]], Fig_I[[2]], Fig_I[[3]], Fig_I[[4]], 
                    align = "hv", labels = LETTERS, font.label = list(size = 16),
                    common.legend = T, legend = "top", nrow = 1, ncol = 4, label.y = 1.05)
fig6.3 <- ggarrange(Fig_N[[1]], Fig_N[[2]], Fig_N[[3]], Fig_N[[4]], 
                    align = "hv", labels = LETTERS, font.label = list(size = 16),
                    common.legend = T, legend = "top", nrow = 1, ncol = 4, label.y = 1.05)
fig6.4 <- ggarrange(Fig[[1]], Fig[[2]], Fig[[3]], Fig[[4]], 
                    align = "hv", labels = LETTERS, font.label = list(size = 16),
                    common.legend = T, legend = "top", nrow = 2, ncol = 2, label.y = 1.05)

fig6.5 <- ggarrange(Fig_I[[6]], Fig_I[[7]], Fig_I[[8]], Fig_I[[9]], Fig_I[[10]], Fig_I[[11]], Fig_I[[12]], Fig_I[[13]], 
                    align = "hv", labels = LETTERS, font.label = list(size = 16),
                    common.legend = T, legend = "top", nrow = 2, ncol = 4, label.y = 1.05)
fig6.6 <- ggarrange(Fig_N[[6]], Fig_N[[7]], Fig_N[[8]], Fig_N[[9]], Fig_N[[10]], Fig_N[[11]], Fig_N[[12]], Fig_N[[13]],
                    align = "hv", labels = LETTERS, font.label = list(size = 16),
                    common.legend = T, legend = "top", nrow = 2, ncol = 4, label.y = 1.05)
fig6.7 <- ggarrange(Fig[[6]], Fig[[7]], Fig[[8]], Fig[[9]], Fig[[10]], Fig[[11]], Fig[[12]],  Fig[[13]],
                    align = "hv", labels = LETTERS, font.label = list(size = 16),
                    common.legend = T, legend = "top", nrow = 4, ncol = 2, label.y = 1.05)

fig6.8 <- ggarrange(Fig_I[[5]], Fig_I[[14]], Fig_I[[15]], Fig_I[[16]], 
                    align = "hv", labels = LETTERS, font.label = list(size = 16),
                    common.legend = T, legend = "top", nrow = 1, ncol = 4, label.y = 1.05)
fig6.9 <- ggarrange(Fig_N[[5]], Fig_N[[14]], Fig_N[[15]], Fig_N[[16]],
                    align = "hv", labels = LETTERS, font.label = list(size = 16),
                    common.legend = T, legend = "top", nrow = 1, ncol = 4, label.y = 1.05)
fig6.10 <- ggarrange(Fig[[5]], Fig[[14]], Fig[[15]], Fig[[16]], 
                    align = "hv", labels = LETTERS, font.label = list(size = 16),
                    common.legend = T, legend = "top", nrow = 2, ncol = 2, label.y = 1.05)

fig_set <- data.frame(name = c("fig6.1", "fig6.2", "fig6.3", "fig6.4", "fig6.5", 
                               "fig6.6", "fig6.7", "fig6.8", "fig6.9", "fig6.10"),
                      width = c(10,15,15,15,15,15,15,15,15,15),
                      height = c(7,3,3,6,6,6,12,3,3,6))
for (i in c(1:10)) {
  ggsave(file = paste(fig_set$name[i], ".png", sep = ""), 
         plot = eval(parse(text = fig_set$name[i])), 
         width = fig_set$width[i], height = fig_set$height[i],  bg = "white")
}


